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Abstract. Different sets of metastable states can be reached in glassy systems below some transition 
temperature depending on initial conditions and details of the dynamics. This is investigated for the 
Sherrington-Kirkpatrick spin glass model with long ranged interactions. In particular, the time dependent 
local field distribution and energy are calculated for zero temperature. This is done for a system quenched to 
zero temperature, slow cooling or simulated annealing, a greedy algorithm and repeated tapping. Results 
are obtained from Monte-Carlo simulations and a Master-Fokker-Planck approach. A comparison with 
replica symmetry broken theory, evaluated in high orders, shows that the energies obtained via dynamics 
are higher than the ground state energy of replica theory. Tapping and simulated annealing yield on the 
other hand results which are very close to the ground state energy. The local field distribution tends to zero 
for small fields. This is in contrast to the Edwards flat measure hypothesis. The distribution of energies 
obtained for different tapping strengths does again not follow the canonical form proposed by Edwards. 

PACS. 05.70.Ln Nonequilibrium and irreversible thermodynamics - 45.70.Cc Static sandpiles; granular 
compaction - 75.10.Nr Spin-glass and other random models - 89.75.-k Complex systems - 02.60.Pn Nu- 
merical optimization 



1 Introduction and Summary 

Complex disordered systems axe ubiquitous, in physics 
and in many other disciplines. Glasses, spin glasses, gran- 
ular media, structure of proteins, neural networks and 
various combinatorial optimization problems have a com- 
plex organization of low energy states in common. Several 
methods have been developed over the time in order to 
deal with the built in disorder. Most widely used is the 
replica method aiming at the evaluation of the free energy, 
entropy or at zero temperature the complexity [1]. Typi- 
cally solutions with spontaneously broken replica symme- 
try show up below some critical value of temperature or 
external noise. 

Alternatively some kind of stochastic dynamics has 
been employed as a tool to investigate such systems [H 
3,4, 5J. For systems with continuous degrees of freedom 
Langevin dynamics may be used. For systems with dis- 
crete degrees of freedom a master equation is appropriate. 
In particular Glauber dynamics is used for Ising-spins. Dy- 
namical processes of this kind can also be used as algo- 
rithm for optimization problems, for example simulated 
annealing [B]. 

The general picture of complex disordered systems is 
associated with a rough landscape of energy or free energy 
with barriers, some of which diverge with the number N 
of particles or elements sent to oo. Within the approach 
via dynamics the limit N — > oo is typically performed 



first. Only thereafter long time scales are eventually in- 
vestigated. This means that the barriers diverging with TV 
can not be overcome and the system might be stuck in 
a certain region of phase space. This means that replica 
theory, not relying on any kind of dynamics, and the ap- 
proach via dynamics might lead to different results. 

This has consequences for instance in using some kind 
of dynamical process for finding solutions of combinato- 
rial optimization problems. Replica theory might tell that 
perfect solutions exist. Dynamics, typically a polynomial 
algorithm, can indicate that these solutions are not found 
in polynomial time. But it provides information about 
suboptimal solutions which can be found in polynomial 
time. An example, learning in a perceptron with binary 
couplings, has been discussed [3]. 

The dynamical behavior of complex disordered sys- 
tems on long time scales is crucially affected by the ex- 
istence of metastable states [3l[il[5lj^[51|^[iT?l[iTl[El] . Sys- 
tems undergoing a discontinuous ergodicity breaking tran- 
sition, e.g. the spherical p-spin-intcraction spin glass, show 
a freezing temperature which is higher than the temper- 
ature where single step replica symmetry breaking sets 
in. This is due to the existence of a large number of meta- 
stable states. In systems with continuous ergodicity break- 
ing transition and full replica symmetry breaking, e.g. 
in the Sherrington-Kirkpatrick model |13j . both tempera- 
tures are identical. Nevertheless the states reached at low 
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temperatures for long, but finite, time may be different to 
those captured by replica theory. 

Dynamics in this context has been developed along es- 
sentially two different lines. Most investigations are based 
on two time correlation and response functions of the re- 
spective basic degrees of freedom [2,3,4,5. Alternatively 
for spin systems and Glauber dynamics the time depen- 
dent local field distribution has been investigated within 
a combination of a Master- and a Fokker-Planck-equation 
[T^lfTSlfTB] . The two methods are in some sense comple- 
mentary. 

There has been recent interest on metastable states 
in the context of granular media. Edwards et.al [17, f 8 
have postulated that the steady state reached in granular 
media after repeated tapping is described by a flat or bi- 
ased average over all metastable states. The bias is in the 
form of a Boltzmann factor with some effective temper- 
ature characterizing the process of tapping. The number 
of metastable states has been computed earlier by Ed- 
wards and others in an annealed approximation [T^ll2Uj . 
This hypothesis was tested with diverse results for vari- 
ous systems and tapping procedures |21U22j . In particular 
Eastham et.al. [TH] have attributed a failure of the Ed- 
wards hypothesis to the special distribution of metastable 
states selected by dynamics. 

Metastable states are commonly defined as local min- 
ima of a free energy functional obtained for instance by 
the TAP approach [TT|[T2] . Investigations of the neighbor- 
hood of those minima or other stationary points yield in- 
teresting results about the structure of low lying states 
in different models of disordered systems. Staying com- 
pletely within the framework of non equilibrium dynamics 
the concept of a free energy does not apply and one has 
to rely on different criteria, for instance on the stability 
of states which are stable with respect to a certain move 
class. Obviously this is restricted to zero temperature and 
a state which is stable with respect to one move class 
might be unstable with respect to a wider move class. 

The present contribution resumes the question what 
kind of metastable states may be reached by various pro- 
cedures. The analysis is based on the temporal behavior of 
the local field distribution. In particular the Sherrington- 
Kirkpatrick model |13j with single spin flip Glauber dy- 
namics is investigated. This is actually a prototype sys- 
tem for a continuous ergodicity breaking transitions or 
for full replica symmetry breaking. A metastable state in 
the present context is a state where each spin points in the 
direction of its local field. This is the field created by the 
external field and the interaction with other spins. Such a 
state is stable with respect to single spin flip dynamics. 

Section [2] contains the definition of the SK- model, the 
local field distribution and Glauber dynamics, which is 
also the basis of the Monte-Carlo simulations presented 
later. The local field distribution obtained from the Ed- 
wards measure [TBl[T51[2Uj is discussed in Section [3] and 
Appendix A. The local field distribution has also been 
computed in high orders of replica symmetry breaking by 
Oppermann et.al [24 . This is discussed for comparison in 
Section [4] The Master- Fokker-Planck approach is the con- 



tent of Section [5] and Appendix B. Results are presented 
and discussed for various assumptions about the drift ve- 
locity in Sections |5.2| and 5.4 Result s fro m the closure 
proposed in 14] are given in Section 5.5 and compared 
with Monte-Carlo simulations in Section |6J 

The following results are obtained: 

The local field distribution at zero temperature and 
late time behaves as P(k) ~ k for k — > + . This holds for 
the following schedules investigated: Quench from a fully 
magnetized or random initial state, greedy algorithm, ran- 
dom or thermal tapping and simulated annealing. This be- 
havior is also found within multiple step replica symmetry 
breaking theory. It contradicts the Edwards flat measure 
hypothesis. 

The distribution of energies found with tapping follows 
a displaced Gaussian, where offset and width depend on 
the tapping strength. According to the Edwards hypoth- 
esis the width should be constant and only the offset is 
supposed to depend on the tapping strength. This is in 
contrast to the present findings. 

Within the framework of a Master-Fokker-Planck ap- 
proach the drift velocity diverges as v(k) ~ for 
k — > + and late time. An approximative closure of the 
resulting hierarchy yields results similar to those obtained 
from Monte-Carlo simulations at zero temperature. 

The ground state energy obtained from replica theory 
is E/N « —0.763. Quenching and greedy algorithm yield 
energies E/N w —0.729 whereas with repeated tapping 
or simulated annealing E/N sa —0.760 is found, which is 
rather close to the ground state energy. The local field dis- 
tribution is also quite close to the one obtained in replica 
theory This indicates that a polynomial algorithm might 
be able to find good, but suboptimal, solutions to a prob- 
lem where finding the best solution is a NP-problem. The 
actual performance certainly depends on on the kind of 
problem and there might be more efficient polynomial al- 
gorithms. 



2 SK-model and Glauber-dynamics 

The energy of the SK-model is given by 

h = - \ Jaw] - hi(ji - w 



where <7j = ±l are Ising spins. The couplings Jij are ran- 
dom variables with 



Jij Jji 



J ij= JZ = N-\ (2) 



In addition to the external field hi a spin Oi on site i feels 
a contribution due to the interaction with other spins. 
Instead of dealing with the resulting local fields, it is con- 
venient to introduce the product of local field and spin 



h = (hi + Y JijVjJVi 



(3) 
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The central quantity of the present investigation is the 
distribution of those fields, sorted according to the direc- 
tion of the spin <7j = <7=±l 



P*(k) 



N 



S a ,a-jS(k - kj) 



(4) 



The bar indicates average over the couplings Jij and the 
brackets denote average over spin configurations. The 
bonds are assumed to be frozen, but the spin configura- 
tions can change in time, resulting in a time dependent 
distribution P a (k\t). In general P a (k) depends on a due 
to the action of external fields or non symmetric initial 
conditions. 

Unless mentioned otherwise stochastic single spin flip 
Glauber dynamics is investigated. At a temperature 
T= 1//3 the flip rate for a spin at site i is 



1 — tanh(/3fci 



(5) 



This is essentially the dynamics used in Monte-Carlo sim- 
ulations and the time unit corresponds to MC-step per 
site. 

In particular at temperature T = the flip rate van- 
ishes for k > and eventually a stationary state with 
ki > for all i is reached. Depending on initial condi- 
tions and cooling schedule various metastable states are 
reached. The distribution of those states depends on initial 
conditions and cooling schedule. 



3 Edwards measure 

Tanaka and Edwards [151120] estimated the number of sin- 
gle spin flip metastable states, i.e. states with ki > 0, in 
annealed approximation. They find for the SK-model a 
displaced Gaussian distribution of energies of metastable 
states 

PMs(e) = \f¥U-^ 2 (6) 



2ir5 



where e = E/N is the energy per site. They estimated 
era 0.5 and 6 = 0.31. 

It has been argued [171115] that a flat or biased aver- 
age over all metastable states applies to tapped granular- 
systems or spin glasses. More precisely, the distribution of 
energies of metastable states obtained with some kind of 
tapping procedure is supposed to behave as 



P E (e) ~P MS (e)e 



(7) 



where the effective temperature fltap depends on the 
strength of tapping. Adopting ([6]), Pfi(e) is again a shifted 
Gaussian with e — > e — S 2 fttap- 

The local field distribution in this approximation, in 
absence of an external field, is a shifted Gaussian trun- 
cated at negative values of k, see [16] and Appendix A, 
Eq. ([501: 



P(fc) = ^P ff (fc) 



-i(fc-«) 2 



erfi-K 



0{k) (8) 
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Fig. 1. Energy e per site (511, local field distribution P(0 + 



and position k of the maximum of P(k) as function of the 
effective inverse temperature fltap calculated from the Edwards 
measure. 




Fig. 2. Local field distributions: a) Edwards measure, Eq. 
b) multi step RSB [21]; c) Master- Fokker-Planck equation; 
d) simplified Master-Fokker-Planck equation; e) Monte Carlo 
simulation at zero temperature; /) Monte Carlo simulation 
with random tapping. 



with 



K= \{P{Q + )+f3ta P }- 



(9) 



This distribution has a discontinuity of size P(0 + ) at k = 0. 
The energy per site is e = — ^[P(0 + ) + «]. 

Eqs. |8| and ^ can be solved numerically for fc = + 
as function of the bias /3 tap . The resulting energy e, discon- 
tinuity of the local field distribution P(0 + ) and peak po- 
sition k are shown in Fig.[T] In Fig. [2] the local field distri- 
bution Q for = 2.5 is plotted together with other results 
discussed later. For the above bias e=— 0.76, k=1.35 and 
P(0 + ) = 0.175. The bias is choosen such that the energy 
e is close to the ground state energy obtained in replica 
theory [24]. 
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4 Multi step replica symmetry breaking 

A linear behavior with slope 0.3 of the local field dis- 
tribution for small fields at zero temperature has been 
found by Sommers and Dupond [23 . Recently Oppermann 
and coworkers |24] have investigated the replica symme- 
try breaking solution for T = in high order. They find 
e= —0.763 ■ • • and again a linear behavior of the local field 
distribution with slope 0.3. The same values have been 
obtained by Pankov [25] ■ The present investigation deals 
with the distribution of the product of local field and spin 
fej, Eq. ([3|, resulting in P{k) — > 0.6 A: for k — > + adopt- 
ing the above value. The complete form of the local field 
distribution is also shown in Fig. [2] 

Even for a bias chosen such that the energy agrees for 
both approaches, P(k) is quite different, indicating that 
the replica calculation and the Edwards measure refer to 
different states. This is actually expected since the replica 
calculation is supposed to present an average over true 
ground states, whereas the Edwards measure is supposed 
to be a biased average over all metastable states. 



5 Master-Fokker-Planck equation 

An equation of motion for the local field distribution has 
been derived in |14j . The following contains a slightly sim- 
plified version. 

A spin flip tji — > —<7i at site i results in changes of the 
modified local fields 



Eq. (12 1 is written as 



and for j ^ i 



kj, > ki 



kj ^ k j 2 J i j ^ % @ ' j • 



(10) 



(11) 



With a flip rate r(k), Eq. (l5j, the resulting time depen- 
dence of the local field distribution is 



d t P a {k, t) = r(-fc)P_ CT (-fc) - r{k)P a {k) + (12) 
^^2\ r ( k i) s °,vi{ s ( k - h + IcnJijOj) - 5(k - kM 



where the first two terms are due to (10 1 and the last 



term is due to (111. Performing the average over the cou- 
plings, neglecting contributions of order -/V -1 , introducing 
a "diffusion constant" (see later) 



d t P a (k,t) = r(-fc)P_ a (-M) -r(fc)P„(M) 
+ 29 fc ^y"dfc' R a<J ,{k,k';t)r{k') 

+ D(t)d 2 k P a (k 7 t). (15) 

The double bracket (•••)) indicates average over P CT (fc, i). 
Introducing a drift velocity v a (k) 



v a {k)P a {k) = -2^2 J dfc'Pw(M'Mfc') 



(16) 



Eq. ([15| combines elements of a Master-equation and a 
Fokker-Planck equation: 



d t P a (k,t) = r(-fc)P_ CT (-M) - r{k)P a {k,t) 
- d k (v a {k,t) - D(t)d*)P a (k,t). 

The two point function ([Til) obeys 



kP tr (k)=Y,[dk' Raa>{k,k') 



and with ( 16 1 the following sum rule is obtained 
(«> = -2(*r). 



(17) 



(18) 



(19) 



Along the same way an equation of motion can be 
derived for the two point function (14 1. The derivation 



is sketched in Appendix B. It involves, however, a three 
point function, and ultimately a hierarchy of equations is 
generated which requires some kind of truncation. 

The probability to flip a spin in field k per unit time 
is r(k). Let r(t) be the probability that a given spin has 
flipped within time t. Then 



d t T{t) = Y J J^r{k)P (J {k,t)=\D{t). 



(20) 



This quantity has actually been used in |16j as a measure 
of time. 



5.1 Initial state 

For any factorizing initial state which is not correlated 
with the couplings Jij , the average over J in (|14j involves 
only the J-dependence of the fields fc,-, Eq. (|3f7 resulting 



D(t) = 2(r{t)} = 2j2j d kr(k)P a {k,t) (13) R aa .{k,k') 



and a two point function 

PfTcr' (^j k ) = 



~N X] ( S °^i S ( k ~ h^iJijCTjS^^Sik' - kj) j (14) 



~ - (d k + <9fc') \ S <r,<n S ( k ~ k i) 5 T',<r j S ( k ' - k ]) 

= -(d k + d k ,)p a {k)P a ,{k>). (21) 

Examples are a fully magnetized initial state with cr; = 1 
for all i, or a state with random spins. In absence of an 
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external field a fully magnetized initial state is actually 
equivalent to any other initial spin configuration within 
the ensemble of couplings. The magnetization plays then 
the role of the overlap with this initial state. 

In the following a fully magnetized initial state is as- 
sumed. The initial value of the local field distribution is 



Especially for t — > oo one obtains 



P+i(M) 



2jt 



P_i(M) = 0. (22) 



The diffusion constant, Eq. (131, is D{Q) — 1 and the drift 
velocity obtained from (16 1 and (21 1 is 



IV (M) - -2«fcr(Q)» 



5.2 Constant and linear drift velocity 



(23) 



Investigating the validity or failure of the Edwards hy- 
pothesis [T5UTT] Eastham et al. [TB] have assumed a con- 
stant drift velocity v(k, t) = c D(t) and have solved Eq. ( 17 1 
numerically for T = and time up to t ~ 10. At t = 10 
they find reasonable over all agreement with MC-simula- 
tions. There appears, however, a slightly smeared out step 
at k = 0. 

The asymptotic behavior at late time can easily be 
evaluated for v{k) — (c—c'k) D(t). Neglecting for a moment 
the non local contributions in (17 1, a Gaussian centered 
at k = n(t) and with width A k {t) results 



Introducing r(t)= L D(t)dt 



the width is 



(l-(l-c')e- 2cV ). 



(24) 



(25) 



(26) 



For T = and late time the part of the local field 
distribution extending to k < is expected to be small. 
This means that D(t), Eq. (fl3j ), is small. For k < the 
dominant contributions to ( |17[ ) are the second term and 
the term involving the second derivative, i.e. 



- P(k, t) + D{t)d 2 k P(k, t)si0 for k < 
which is solved by 

P(M) = a(t)e k /v^, 



(27) 



(28) 



and with (13 1 D(t) — 2a 2 (t). Similar arguments lead to 
P(k,t) = a(t)(2-e- k /V°W} 



(29) 



for small k > 0. The factor a(t) is obtained by matching 
with the solution (24) at k ~ 0, i.e. 



a(t) 



1 



2 v /2 7 rZ\(t) 



-K 2 (t)/2A 2 (t) 



(30) 



c 
c 



A 



1 



D 



r j- 



8tt 



-c 2 /c' 



(31) 



This solution resembles the in some sense the distribu- 
tion obtained from the Edwards measure ^ with (3 = 2c 
and d = 1. The step at k = is, however, smeared out 
and even at late time there is a tail of P(k) extending to 
k < indicating that the above assumption regarding the 
drift velocity is not appropriate for metastable states with 
P(k < 0)=0. 

For d — > and t — * oo the step at k — vanishes, but 
K and with it the energy per site grow without limit. The 
qualitative agreement with the MC-data in [T5] appears to 
be a consequence of the particular choice of c and an ap- 
propriate finite t in this investigation. It has been argued 
that the results at finite time are actually more appropri- 
ate for the asymptotic behavior of a system of finite size. 
I shall come back to this point later. 



5.3 Qualitative discussion of the drift term at T=0 

At zero temperature and for late time metastable states 
are reached, this means that all local fields are expected 
to be non negative, i.e. P(k < 0, t — > oo) — ■> 0. Assume 
that for small k and late time 



P a (k,t) -> ck a O(k). 



(32) 



This gives rise to a contribution ~ k a "in the diffusion 



term of Master-Fokker-Planck equation ( 17 1, which has to 



be compensated by the drift term. This results in 

v a (k,t) -> aD^r 1 . (33) 

The evaluation of the actual value of a requires to solve 
the equations for the n-point functions of higher order. 

5.4 Closure I: Disregarding dynamical correlations 

The two point function ( |14[ ) involves an average over the 
couplings Jij. The couplings are contained in the local 
fields ki, ([3). In addition correlations between the cou- 
plings and the actual states reached at finite t build up in 
the course of time. Neglecting those Eq. (21 1 holds for all 
times. The resulting drift velocity is 



v a (k,t) = D(t)d k lni y P a (k)j-2(r / (t)j (34) 
and the equation of motion (|17|) becomes 



d t P a {k,t) w r{-k)P_ a {-k,t)-r{k)P a {k,t) 

+ 2{r'(t)}d k P a (k,t). (35) 

The qualitative behavior for T = is easily discussed. 
For small k the local field distribution is 



P a (k,t)^P a (0,t) + kPU0 ± ,t) 



(36) 
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with P+(0,ty f P-{0,t) and P;(0+,i) ^ f£(G~,i). 
fc -> Eq. (pi yields 



For 



9 t ^P CT (0,t) (37) 

(T 

= -2 £ ^(0, f) X! (° + . *) - ^ > *) ' 



which is solved by ^2 -fo-(O)^) ~~ * f° r t ^ oo with finite 
i*(0+,t) > 0. The drift velocity obeys for fc > 



D(t) 
fc 



(38) 



and the exponent a introduced in the previous subsection 
is a = l. 

This is supported by the numerical integration of 
Eq. (35 1. The complete local field distribution for t^oo 
is also shown in Fig. [2] The slope of P(k) for fc — * is 
now 0.9. This increased value is in accordance with the 
maximum of P(k) shifted towards smaller values of fc and 
a higher value of the energy e. This behavior is also found 
in the investigations reported later. 

Already this rather simple closure shows that the as- 
sumption of a constant or linear drift velocity, |16j and 
Sect. |5.2[ is not appropriate. It also shows that the Ed- 
wards measure |19j does not apply to the sit uati on cap- 
tured by the Master-Fokkcr-Planck equation (17). It has 



to be stressed again, that this equation describes a rapid 
quench from a fully magnetized or high temperature state 
to T = 0, and that the limit N— >oo is performed first. In 
terms of combinatorial optimization problems this corre- 
sponds to a greedy search and therefore to a polynomial 
algorithm ~ iV 2 . An average over all metastable states, on 
the other hand, would require an exponential effort. 



5.5 Closure II: Approximative treatment of dynamical 
correlations 

An improved theory taking into account dynamical corre- 
lations has been proposed by the author [T3]. It is based 
on a factorization of the three point function entering the 
equation of motion of the two point function (14 1. It in- 



volves a modified Kirkwood superposition approximation 
known from the theory of real gasses [26 . A slightly sim- 
plified version is outlined in Appendix B. The Master- 
Fokker-Planck equation (17) and the equation for the drift 



velocity ( 55 1 can be integrated numerically with the ini- 
tial conditions discussed in Sect. 15. II For technical reasons 
the calculations are performed for a small finite temper- 
ature, typically T = 0.01, much smaller than the freezing 
temperature T c = 1 . 

Fig. [3] shows the resulting local field distributions for 
spin up and down, respectively, for times 0, 1, 5 and 30. 
The result from the Monte-Carlo simulation discussed in 
the next section for t = 30 are also shown. The complete 
local field distribution P(k) — P + (k) + fL(fc) is compared 
with other results in Fig. [2] This shows that for long time 
P(fc, i) — > for fc — > in qualitative agreement with 




Fig. 3. Local field distribution P a (k,t) for t = 0, 1, 5 and 30 
obtained from the Master-Fokker-Planck approach and from 
Monte-Carlo simulations for t — 30. 
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Fig. 4. Drift velocity v a {k,t) for t = 0, 1, 5 and 30 obtained 
from the Master-Fokker-Planck approach and from Monte- 
Carlo simulations for £ = 30. 



replica theory and in contrast to the finite step obtained 
from Edwards hypothesis. The value of the exponent a 
introduced in Sect. [53] is consistent with a—1. 

The agreement between Monte-Carlo simulation and 
the Master-Fokker-Planck approach with the present clo- 
sure is quite satisfactory. This indicates that the present 
closure captures the essential mechanisms: the effect of 
flipping a spi n in its own local field, described by the first 
two terms in (17 1, and the small effects on the local fields 



of all the other spins, described by the drift- and diffusion 
terms of the Master-Fokker-Planck equation (17 1. 



The drift velocity for spin up and down for the same 
times as above is plotted in Fig.[4]together with the results 
of the simulation at the latest time. The formation of the 
1/fc-divergence of v a (k,t) for long time and fc — > is seen 
in the Master-Fokker-Planck approach as well as in the 
Monte-Carlo data. 

Energy per site e(i), magnetization m(t), diffusion con- 
stant D(t) and spin flips per site r(t), see Eq. (20 1, are 
plotted as functions of time in Fig. [5] Again Master-Fokker- 
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Fig. 5. Energy per site e(t), magnetization m(t), diffusion con- 
stant D(t), flips per site r(t) and fraction To(t) of spins fliped at 
least once obtained from the Master-Fokker-Planck approach 
and from zero temperature Monte-Carlo simulations. The time 
unit corresponds to Monte-Carlo steps per site. 



Planck and Monte-Carlo results are compared. Asymp- 
totic values for t — > oo are listed in Table [2j They are 
discussed later. 



6 Monte-Carlo simulations 

In order to test the results of the previous section and in 
order to investigate alternative optimization methods var- 
ious Monte-Carlo simulations have been performed. The 
general procedure is standard. A site i is selected at ran- 
dom and flipped with a probability r(k j), Eq. |5]). The 
local fields are updated according to (111 



(39) 



The local field distribution, Eq. Q , two point function, 
Eq. (14 1 and drift velocity, Eq. (16 1, are also updated ac- 



cordingly. 

Alternatively a greedy algorithm is investigated, where 
in each step the spin cr, withe most negative value of ki is 
selected and nipped, until ki > for all i. 

For the SK-model in the limit N ~ + oo only the first 
and second moment of the distributions of couplings Q 
is relevant. The simulations are therefore performed with 



Jij =±l/y/N for i ^= j and J a = 0. This speeds up the 
computation considerably and the local fields are integer 
multiples of 1 / y/N. 

The simulations are performed on samples with up to 
N = 22500 sites. The results are typically averaged over 
1000 runs with 200 different sets of random couplings Jij. 





F 


F' 


CMC 


-0.729 


0.157 


^greedy 


-0.730 


0.182 


rriMC 


0.074 


0.29 


TTlgrccdy 


0.091 


0.53 



Table 1. Coefficient of fitting energy per site —e(N) and mag 
netization m(N) according to (40 1 with a = 0.2. 



6.1 Zero temperature simulation and greedy algorithm 

Results of the zero temperature Monte-Carlo simulations 
and from the greedy algorithm are shown in Figs. [2] to [5] 
and in Table El 

Fig.[5]compares the time dependence of various quanti- 
ties obtained from zero temperature Monte Carlo simula- 
tions on samples of different size N and from the Master- 
Fokker-Planck approach. Regarding the energy, magneti- 
zation and diffusion constant there are some discrepancies 
at intermediate times, the behavior at early and late time 
is, however, reasonably well reproduced. 

Energy and magnetization can be fitted according to 



F(N) « F + F' N~ 



(40) 



In [27] a = 0.33 was proposed. An improved fit to the 
present data is obtained with a — 0.2. This yields the val- 
ues shown in Tab.[T| The asymptotic values of the energy 
are slightly lower than those obtained in [37] ■ This is due 
to the different choice of a and the fact that larger values 
of N are used in the present investigation. 

Starting with a fully magnetized initial state the mag- 
netization at T = remains finite. This is also true for 
the greedy algorithm. The values obtained by the simula- 
tions are in reasonable agreement with the result from the 
Master-Fokker-Planck approach. 

The local field distribution P(k) obtained from the 
greedy algorithm and from zero temperature Monte-Carlo 
simulations are almost identical. For finite N there is a 
small step at k — vanishing as 



P(k = + ,N) « 1.4 AT 



-0.46 



(41) 



which is in reasonable agreement with the result obtained 
in [27]. 

In the limit N — > oo P±(k) ~ k for k — * is found 
in all three cases. This is in contrast to the step resulting 
from the Edwards hypothesis. 

The distribution of the energy per site Pe(c), see ([6]), 
is fitted well by a shifted Gaussian, with e w —0.70 and 
5 1=3 0.62 for TV = 10000. Similar values are found for other 
values of N. 

A significant discrepancy between the results of the 
Master-Fokker-Planck approach and Monte Carlo simula- 
tions shows up for the number of spin flips measured as 
r(t). A given spin might actually flip more than once. Let 
To(t) be the fraction of spins which has flipped at least 
once. Then the fraction of spins which have not flipped at 
all is 1 — To(t). This quantity is also shown in Fig. [5] The 
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difference between r(t) and To(t) is a measure of multi- 
ple spin flips. The results obtained for samples of different 
size follow universal curves up to some iV-dependent time 
i(N). A tit covering the range t > 10 yields 

t(N) ps 0.22 TV - 7 (42) 

and for r(t,N) and To(t,N), respectively, the asymptotic 
values 

t(oo, N) ps 0.26 + 0.077 ln(iV) 

r(t,oo) fa 0.43 + 0.11 ln(i) 
t (oo, N) PS 0.66 - 0.19 N-"- 25 
T (t, oo) Ps 0.66 - 0.32 1~ 36 (43) 

are obtained. This means that r(t,n) diverges logarithmi- 
cally for t — > oo and N — > oo whereas To(t,N) remains 
finite in this limit. A divergence of r(t, N) was also found 
by Eastham et.al. [16] on the basis of a constant drift ve- 
locity. A fit corresponding to Eq. (43 1 would, however, 
yield rather different values. The Master-Fokker-Planck 
approach yields a finite asymptotic value for r(i^oo) 
which is not too far from the asymptotic value of t (oo, oo) . 



6.2 Tapping Dynamics 

An investigation of tapping dynamics is of interest not 
only in the context of the Edwards hypothesis [17] • Tap- 
ping or thermal cycling might also be used for combina- 
torial optimization problems |28j . Tapping might be done 
either by heating the system periodically to some tem- 
perature T tap or by flipping randomly selected spins with 
probability p. 

The minimal energy obtained within n tap cycles of ran- 
dom tapping is fitted to 

£ta P (n t ap,P, N) = e 0: ta P (p, N) - e' tap (p, N) ra^p (44) 

with a' = 0.5. There is no significant iV-dependence in 
eo t tap{p, N) and the optimal value eo,tap(P =0-1) ~ —0.759 
is obtained for p — 0.1. 

Thermal tapping yields similar results with an opti- 
mal tapping temperature T tap ps 0.7. Tapping is obviously 
quite effective in finding states with low energy. 

For comparison simulated annealing has also been 
tested. A schedule T(t) = (1 - t/t) T with T = 1.5 and 
i = 300 • • • 10000 has been used. For the slowest cooling 
schedule e ps —0.758 has been found. 

The local field distribution obtained with the tapping 
procedure is shown in Fig. [2] The agreement with the lo- 
cal field distribution obtained by multi step RSB is quite 
good. Again there is a small step of P(k,N) at k = 0, 
which vanishes for N — > oo similar to Eq. (41 ) with a re- 
duced pref actor ^0.5. Even at finite N P{0 + ) is much 
smaller than P{0 + ) ~0.18 resulting from the Edwards flat 
measure hypothesis. 

As a second test of the validity of the Edwards flat 
measure hypothesis the distribution of energies of meta- 
stable states Pms(e) for different tapping strengths are 
compared [22] . 



Method 


— e 


m 


Multi step RSB 


0.763 




Eswards measure /3 = 2.5 


0.760 




Master-Fokker-Planck eq. 


0.688 


0.086 


Monte Carlo T = 


0.729 


0.074 


Greedy 


0.730 


0.091 


Random tapping, p=W% 


0.760 




Thermal tapping, Tt ap = 0.7 


0.759 




Simulated annealing 


0.759 





Table 2. Energy per site e and magnetization m obtained from 
the various methods discussed in the text. For the simulations 
extrapolated values for TV — * oo are shown. 




0.745 



0.750 



0.755 



0.760 



Fig. 6. Integrated distribution of energies #(e) for different 
tapping strengths p using random tapping wit TV = 2500 and 
n = 1000. 



Counting how often an energy Ei — Nei < Ne is found 
within L Monte Carlo runs, the integrated distribution of 
energies of metastable states 



1 f°° 
${e) = jY j G{e-e l )= / de' P MS 



(e') (45) 



results. This quantity is shown in Fig.[6]obtained from ran- 
dom tapping with strength p on systems of size ^ = 2500. 

Adopting a Gaussian distribution of energies of meta- 
stable states, Eq. the integrated distribution is 



#(e) 



(46) 



Combinig this with Eq. ^ , which is a consequence of the 
Edwards hypothesis, the only effect of tapping should be 
a shift of the peak energy e, whereas S should be constant. 
The data in Tab. show that this is not the case. 

The effective temperature 1 / (3tap in Eq. (|7| is supposed 
to depend on the tapping strength p. More precisely 



S 2 Pt ap = const. 



(47) 



is expected. The existence of a finite optimal tapping 
strength implies that this dependence can not be mono- 
tonous. For thermal tapping similar distributions are found. 
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p % 


— c 


8 




— e 


5 


5.0 


0.7509 


0.263 


0.4 


0.7458 


0.227 


7.5 


0.7537 


0.207 


0.5 


0.7495 


0.174 


10.0 


0.7548 


0.163 


0.6 


0.7517 


0.151 


12.5 


0.7544 


0.157 


0.7 


0.7532 


0.138 


15.0 


0.7526 


0.172 


0.8 


0.7536 


0.127 


17.5 


0.7504 


0.188 


0.9 


0.7532 


0.128 


20.0 


0.7468 


0.189 


1.0 


0.7512 


0.134 



with k= 1{P(0 + ) + /?}. The energy e per site is 



Table 3. e and S obtained from a fit of <P(/e) to Eq. (461 for 
random and thermal tapping, respectively. 



Again 6 is not constant and Eq. (47 1 is not fulfilled assum- 
ing fitap = l/T tap . 



Appendix A: Flat average over metastable 
states 

Performing a flat or biased average over all metastable 
states, the local field distribution, say at site o is 



P(k) = Z- l Tx ao []Tr CT . S(k - k )e^o/2 

J 

x Yl 0(h)e^/ 2 O(k) (48) 

ifro) 

with k a and ki given in The metastable states are 
weighted with a Boltzmann factor er^ E where 
E = — | ki. The effective temperature 1//3 is assumed 
to characterize the tapping procedure. It might as well be 
viewed as Lagrange multiplier selecting the total energy 
of the metastable states under consideration. 

The local field at a site i ^ is writen as ki = k' i + Ki 
with Ki = o J iOi . k'i is the local field of the system without 
the spin a . For the SK-model in the limit N — ► oo 

J POO 

Y[ G{k[ + Ki) « \[ / dkiP(h) (49) 
« II {i + O(-k 1 )[k 1 P(0 + )-^P'(0+))}. 

Using the Fourier representation of the 8- function and per- 
forming the average over Ki = ±-k 



N 



P(k) = z- 1 



dk 
2?r ' 



,(fc+/3/2)fc 



{l + ^(fc 2 -fcP(0+)-ip'(0+))} O(k) 

9(k) 



e 2 



i[fc-i{P(0+)+/3}] 2 



- l + erf(^{P(0+)+/3}) 
p — \ (k—K) 

- l + erf(^«) 



(50) 



\ JdkkP(k) = -±(p(0+) + K 



(51) 



Appendix B: Equation of motion for the two 
point function and drift velocity 



Let R^g., (k, k') be the contribution of sites i and j to ( 14 ) 
Flipping cTj — > — (Ji yields 



AR^,(k,k') = - r {-k)R% a ,{-k,k')-r{k)Rt.{Kk') 
- 2r{-k)P- <r {-k)d k ,P a ,{k'). (52) 

The first two terms are due to — fcj — ► ki. The last term 
has its origin in kj — > fcj — 2uiJijGj. Actually this term 



-./' 



with Jy = 0. 



contains (8- a _ ai S(k + ki)S a ' . aj S(k' — kj)) 
This expression factorizes, however, for the SK-model. A 
corresponding contribution arises from flipping spin cry. 

Flipping a spin on site I ^ i,j yields diffusion terms 
and drift terms in analogy to (151. The latter actually 
involve three point functions. In~^14] a closure has been 
proposed on the basis of a generalized Kirkwood superpo- 
sition approximation. A slight simplification yields 



£ ARZ,(k,k') a - d k (v a {k) - D{t)d k ) 

+[k^k']\RZ,(k,k'), (53) 



Collecting the above contributions 
d t R aryl (k,k') « -r{-k)R^(-k,k') - r{k)R aa ,{k,k') 
- dkfv^-Dd^R^i^k') 
-2r(-k)P_ a (-k)d k/ P a ,(k>) 



a', k^k' 



(54) 



This yields the equation of motion for the drift velocity 
( fl6| . There are different types of contributions. The first 
two lines of ( |54| result in corresponding expressions for 
v a (k)P a (k). The contribution of the third line is easily 
taken into account, because it factorizes. The terms anal- 
ogous to the first two lines acting on a' and k' can not 
be expressed in terms of P and v only. Taking into ac- 
count the identities (151 and ( 18 1 they are approximated 
by adding contributions ~ v a (kJP a (k) and ~ kP a (k). The 
resulting equation of motion for v a (fc) is 

d t v a (k) = -r(-k) (y a (k) + «_,(-*) + 4 (r'») 



v a {k) - Dd k - 2D\d k \n{P a {k))^d k v a {k) 

+ 4(rf) \d k \n(P a (k)) \ +Av a (k) + Bk. (55) 

For k — > oo the local field distribution is expected to follow 
P a {k > 1) - e~5 fe2 . This is the case if (EH) holds not 
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only for the initial state, but for all time in the limit k or 
kl — » oo. The coefficients A and B are then determined 



such that this asymptotic behavior and the sum rule (19 1 
is fulfilled. 



I like to thank Mike Moore for discussions and Paul Eastham 
for helpful correspondence, especially on the Edwards hypoth- 
esis and the resulting local field distribution. 
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